function [J, C] = computeQuadrotorMatrices(phi, theta, psi, phi_dot, theta_dot, psi_dot, I11, I22, I33)
    % Inputs:
    %   phi, theta, psi       - Euler angles (roll, pitch, yaw) in radians
    %   phi_dot, theta_dot, psi_dot - Angular velocities (rad/s)
    %   I11, I22, I33         - Inertia moments (kg·m^2)
    %
    % Outputs:
    %   J                     - Inertia matrix
    %   C                     - Coriolis matrix

    % Calculate inertia matrix J
    J11 = I11;
    J13 = -I11 * sin(theta);
    J22 = I22 * cos(phi)^2 + I33 * sin(phi)^2;
    J23 = (I22 - I33) * cos(phi) * sin(phi) * cos(theta);
    J33 = I11 * sin(theta)^2 + I22 * sin(phi)^2 * cos(theta)^2 + I33 * cos(phi)^2 * cos(theta)^2;

    J = [J11, 0,    J13;
         0,    J22, J23;
         J13, J23, J33];

    % Calculate Coriolis matrix C
    c11 = 0;

    c12 = (I22 - I33) * (theta_dot * cos(phi) * sin(phi) + psi_dot * sin(2 * phi) * cos(theta) - psi_dot * cos(phi)^2 * cos(theta)) ...
          - I11 * psi_dot * cos(theta);
    c13 = (I33 - I22) * psi_dot * cos(phi) * sin(phi) * cos(theta)^2;

    c21 = (I33 - I22) * (theta_dot * cos(phi) * sin(phi) + psi_dot * sin(2 * phi) * cos(theta) - psi_dot * cos(phi)^2 * cos(theta)) ...
          + I11 * psi_dot * cos(theta);
    c22 = (I33 - I22) * phi_dot * cos(phi) * sin(phi);
    c23 = -I11 * psi_dot * sin(theta) * cos(theta) ...
          + I22 * psi_dot * sin(phi)^2 * cos(theta) * sin(theta) ...
          + I33 * psi_dot * cos(phi)^2 * cos(theta) * sin(theta);

    c31 = (I22 - I33) * psi_dot * cos(theta)^2 * sin(phi) * cos(phi) ...
          - I11 * theta_dot * cos(theta);
    c32 = (I33 - I22) * (theta_dot * cos(phi) * sin(phi) * sin(theta) + phi_dot * sin(2 * phi) * cos(theta) - phi_dot * cos(phi)^2 * cos(theta)) ...
          + I11 * psi_dot * sin(theta) * cos(theta) ...
          - I22 * psi_dot * sin(phi)^2 * sin(theta) * cos(theta) ...
          - I33 * psi_dot * cos(phi)^2 * sin(theta) * cos(theta);
    c33 = (I22 - I33) * phi_dot * cos(phi) * sin(phi) * cos(theta)^2 ...
          - I22 * theta_dot * sin(phi)^2 * cos(theta) * sin(theta) ...
          - I33 * theta_dot * cos(phi)^2 * cos(theta) * sin(theta) ...
          + I11 * theta_dot * cos(theta) * sin(theta);

    C = [c11, c12, c13;
         c21, c22, c23;
         c31, c32, c33];
end
